library(haven)
library(miceadds)
library(sandwich)
library(dplyr)
library(ggplot2)
library(tidyr)
library(ggpubr)
library(scales)
library(interplot)
library(dotwhisker)


# Open the data set called "data.dta" and save it as a data frame called "deservingness"
#deservingness <- read_dta("...data.dta")

# Recoding 
deservingness$treat_corona <- as.factor(deservingness$treat_corona)

# Descriptive statistics (Table 3)
treatmentgroup <- deservingness %>%
  filter(treat_corona==1)

controlgroup <- deservingness %>%
  filter(treat_corona==0)

# Column 1 ("Total") of Table 3 
mean(deservingness$kon, na.rm = TRUE)
sd(deservingness$kon, na.rm = TRUE)

mean(deservingness$alder, na.rm=TRUE)
sd(deservingness$alder, na.rm=TRUE)

mean(deservingness$erfaring, na.rm=TRUE)
sd(deservingness$erfaring, na.rm=TRUE)

mean(deservingness$brugerorientering, na.rm=TRUE)
sd(deservingness$brugerorientering, na.rm=TRUE)

mean(deservingness$gruppe, na.rm=TRUE)
sd(deservingness$gruppe, na.rm=TRUE)

# Column 2 ("Control") of Table 3
mean(controlgroup$kon, na.rm = TRUE)
sd(controlgroup$kon, na.rm = TRUE)

mean(controlgroup$alder, na.rm=TRUE)
sd(controlgroup$alder, na.rm=TRUE)

mean(controlgroup$erfaring, na.rm=TRUE)
sd(controlgroup$erfaring, na.rm=TRUE)

mean(controlgroup$brugerorientering, na.rm=TRUE)
sd(controlgroup$brugerorientering, na.rm=TRUE)

mean(controlgroup$gruppe, na.rm=TRUE)
sd(controlgroup$gruppe, na.rm=TRUE)

# Column 3 ("Treatment") of Table 3
mean(treatmentgroup$kon, na.rm = TRUE)
sd(treatmentgroup$kon, na.rm = TRUE)

mean(treatmentgroup$alder, na.rm=TRUE)
sd(treatmentgroup$alder, na.rm=TRUE)

mean(treatmentgroup$erfaring, na.rm=TRUE)
sd(treatmentgroup$erfaring, na.rm=TRUE)

mean(treatmentgroup$brugerorientering, na.rm=TRUE)
sd(treatmentgroup$brugerorientering, na.rm=TRUE)

mean(treatmentgroup$gruppe, na.rm=TRUE)
sd(treatmentgroup$gruppe, na.rm=TRUE)

# Test of differences between groups
summary(aov(kon ~ treat_corona, data=deservingness))
summary(aov(alder ~ treat_corona, data=deservingness))
summary(aov(erfaring ~ treat_corona, data=deservingness))
summary(aov(brugerorientering ~ treat_corona, data=deservingness))
summary(aov(gruppe ~ treat_corona, data=deservingness))

# Table 4
m1 <- lm.cluster(manipulationstjek ~ treat_corona, cluster = "kommune", data=deservingness)
summary(m1)

m2 <- lm.cluster(corona_sanktion ~ treat_corona, cluster = "kommune", data=deservingness)
summary(m2)

m3 <- lm.cluster(corona_hjaelp ~ treat_corona, cluster = "kommune", data=deservingness)
summary(m3)

# Figure 1
long_deservingness <- deservingness %>% 
  pivot_longer(c(manipulationstjek, corona_sanktion, corona_hjaelp), names_to ="DV", values_to = "response")

# Helper function 

summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  library(plyr)
  
  # New version of length which can handle NA's: if na.rm==T, don't count them
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
  }
  
  # This does the summary. For each group's data frame, return a vector with
  # N, mean, and sd
  datac <- ddply(data, groupvars, .drop=.drop,
                 .fun = function(xx, col) {
                   c(N    = length2(xx[[col]], na.rm=na.rm),
                     mean = mean   (xx[[col]], na.rm=na.rm),
                     sd   = sd     (xx[[col]], na.rm=na.rm)
                   )
                 },
                 measurevar
  )
  
  # Rename the "mean" column    
  datac <- rename(datac, c("mean" = measurevar))
  
  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  
  # Confidence interval multiplier for standard error
  # Calculate t-statistic for confidence interval: 
  # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  
  return(datac)
}

m1 <- summarySE(long_deservingness, measurevar="response", groupvars=c("treat_corona", "DV"), na.rm=T) 

dv.labs <- c("Help client", "Recommend sanction", "Responsibility") 
names(dv.labs) <- c("corona_hjaelp", "corona_sanktion", "manipulationstjek")


p <- ggplot(data=subset(m1, !is.na(treat_corona)), aes(x=treat_corona, y=response, fill=treat_corona)) + geom_bar(stat="identity")+ 
  geom_errorbar(aes(ymin=response-ci, ymax=response+ci), width=.05, position=position_dodge(.9)) + theme_bw() +
  labs(x="", y = "Mean response (0-10)") + 
  scale_fill_grey(start=0.6, end=0.3, name="", labels=c("Low controllability", "High controllability")) +
  scale_x_discrete(labels=c("0" = "", "1" = "")) + ylim(0,10) +
  facet_wrap(vars(DV), ncol=3, labeller=labeller(DV=dv.labs)) + theme(legend.position="bottom")

#pdf("fig1.pdf")
#p
#dev.off()


# Figure 2

#pdf("Fig2.pdf")
#par(mfrow=c(1,2))

hist_sanktion <- hist(deservingness$corona_sanktion, plot = FALSE)  # Store output of hist function
hist_sanktion$density <- hist_sanktion$counts /    # Compute density values
  sum(hist_sanktion$counts) * 100
plot(hist_sanktion, freq = FALSE, xlab="Recommend sanction", ylab="Percentage", main="", ylim=c(0,25))  # Plot histogram with percentages

hist_help <- hist(deservingness$corona_hjaelp, plot = FALSE)    # Store output of hist function
hist_help$density <- hist_help$counts /    # Compute density values
  sum(hist_help$counts) * 100
plot(hist_help, freq = FALSE, xlab="Helping client", ylab="Percentage", main="", ylim=c(0,25)) # Plot histogram with percentages

#dev.off()


# Figure 3
kommune_filter_sanktion <- deservingness %>%
  filter(kommune %in% c('2','4','5','10','11','12','14','20','21','24','25','26','27','28','30'))

ggplot(kommune_filter_sanktion, aes(corona_sanktion)) + geom_histogram(aes(y=after_stat(width*density))) + 
  facet_wrap(~kommune, ncol=3) + labs(x="Recommend sanction", y="Percentage") +
  scale_y_continuous(labels = percent_format()) + scale_x_continuous(breaks = c(0,1,2,3,4,5,6,7,8,9,10)) + theme_bw()

# Figure 4
kommune_filter_help <- deservingness %>%
  filter(kommune %in% c('2','4','5','10','11','12','14','20','21','24','25','26','27','28','30'))

ggplot(kommune_filter_help, aes(corona_sanktion)) + geom_histogram(aes(y=after_stat(width*density))) + 
  facet_wrap(~kommune, ncol=3) + labs(x="Helping client", y="Percentage") +
  scale_y_continuous(labels = percent_format()) + scale_x_continuous(breaks = c(0,1,2,3,4,5,6,7,8,9,10)) + theme_bw()


# Figure 6
deservingness$treat_corona.f <- as.factor(deservingness$treat_corona)

#pdf("Fig6.pdf")
#par(mfrow=c(1,2))

mod_hjaelp <- lm(corona_hjaelp ~ treat_corona.f*alder, data=deservingness)
summary(mod_hjaelp)

interplot(m=mod_hjaelp, var1="treat_corona.f", var2 ="alder", facet_labs=(""), hist=T, ci=.90) + theme_bw() + geom_hline(yintercept=0, linetype="dashed") +
  xlab("Age (years)") + ylab("Marginal effect of treatment on willingness to help client") 

mod_hjaelp2 <- lm(corona_hjaelp ~ treat_corona.f*erfaring, data=deservingness)

summary(mod_hjaelp2)

interplot(m=mod_hjaelp2, var1="treat_corona.f", var2 ="erfaring", facet_labs=(""), hist=T, ci=.90) + theme_bw() + geom_hline(yintercept=0, linetype="dashed") +
  xlab("Tenure (years)") + ylab("Marginal effect of treatment on willingness to help client") 

#dev.off()

# Figure 5

# Calculating point estimates and standard errors via regression 
mod_age_hjaelp <- lm(corona_hjaelp~alder + treat_corona, data=deservingness)
mod_age_treat_hjaelp <- lm(corona_hjaelp~alder*treat_corona, data=deservingness)
mod_female_hjaelp <- lm(corona_hjaelp~kon + treat_corona, data=deservingness)
mod_female_treat_hjaelp <- lm(corona_hjaelp~kon*treat_corona, data=deservingness)
mod_tenure_hjaelp <- lm(corona_hjaelp~erfaring + treat_corona + alder + kon, data=deservingness)
mod_tenure_treat_hjaelp <- lm(corona_hjaelp~erfaring*treat_corona + alder + kon, data=deservingness)
mod_user_hjaelp <- lm(corona_hjaelp~brugerorientering + treat_corona + alder + kon + erfaring, data=deservingness)
mod_user_treat_hjaelp <- lm(corona_hjaelp~brugerorientering*treat_corona + alder + kon + erfaring, data=deservingness)

mod_age_sanktion <- lm(corona_sanktion~alder + treat_corona, data=deservingness)
mod_age_treat_sanktion <- lm(corona_sanktion~alder*treat_corona, data=deservingness)
mod_female_sanktion <- lm(corona_sanktion~kon + treat_corona, data=deservingness)
mod_female_treat_sanktion <- lm(corona_sanktion~kon*treat_corona, data=deservingness)
mod_tenure_sanktion <- lm(corona_sanktion~erfaring + treat_corona + alder + kon, data=deservingness)
mod_tenure_treat_sanktion <- lm(corona_sanktion~erfaring*treat_corona + alder + kon, data=deservingness)
mod_user_sanktion <- lm(corona_sanktion~brugerorientering + treat_corona + alder + kon + erfaring, data=deservingness)
mod_user_treat_sanktion <- lm(corona_sanktion~brugerorientering*treat_corona + alder + kon + erfaring, data=deservingness)

# Point estimates
coef.matrix.1 <- matrix(c(0.02, -0.006,
                          0.02, -0.002,
                          -0.03, 0.01,
                          0.06, -0.03), nr=1)

# Standard errors
se.matrix.1 <- matrix(c(0.006, 0.009,
                        0.012, 0.017,
                        0.011, 0.016,
                        0.02, 0.028), nr=1)


# Variable names
varnames <- c("")

model_names <- c("Age", "Age*treatment", "Tenure", "Tenure*treatment")


submodel_names <- c("Help-giving", "Sanctioning")

results_df_1 <- data.frame(term=rep(varnames[1], times=2),
                           estimate=as.vector(coef.matrix.1),
                           std.error=as.vector(se.matrix.1),
                           model=as.factor(rep(model_names, each=2)),
                           submodel=rep(rep(submodel_names, each = 1), times = 4))

explorative_plot_1 <- small_multiple(results_df_1) + 
  scale_x_discrete(limits=model_names) +
  theme_bw() + ylab("Marginal effects") +
  geom_hline(yintercept = 0, colour = "black", linetype = "dashed") +
  theme(axis.text.x = element_text(size=10),
        axis.text.y = element_text(size=10),
        axis.title.y = element_text(size=12),
        legend.position=c(.99, .99), legend.justification=c(1, 1),
        legend.title = element_text(size=12),
        legend.background = element_rect(color="black"),
        legend.key.size = unit(12, "pt")) +
  scale_colour_grey(start=0.2, end=0.8, name="Dependent variable") +
  ggtitle("")

explorative_plot_1

# Point estimates
coef.matrix.2 <- matrix(c(0.005, -0.08,
                          0.59, -0.27,
                          1.18, -1.59,
                          -0.67, 1.99), nr=1)

# Standard errors
se.matrix.2 <- matrix(c(0.17, 0.24,
                        0.35, 0.49,
                        0.42, 0.60,
                        0.85, 1.19), nr=1)

# Variable names
varnames <- c("")

model_names <- c("Female", "Female*treatment", "User orientation", "User orientation*treatment")


submodel_names <- c("Help-giving", "Sanctioning")

results_df_2 <- data.frame(term=rep(varnames[1], times=2),
                           estimate=as.vector(coef.matrix.2),
                           std.error=as.vector(se.matrix.2),
                           model=as.factor(rep(model_names, each=2)),
                           submodel=rep(rep(submodel_names, each = 1), times = 4))


explorative_plot_2 <- small_multiple(results_df_2) + 
  scale_x_discrete(limits=model_names) +
  theme_bw() + ylab("Marginal effects") +
  geom_hline(yintercept = 0, colour = "black", linetype = "dashed") +
  theme(axis.text.x = element_text(size=10),
        axis.text.y = element_text(size=10),
        axis.title.y = element_text(size=12),
        legend.position=c(.99, .99), legend.justification=c(1, 1),
        legend.title = element_text(size=12),
        legend.background = element_rect(color="black"),
        legend.key.size = unit(12, "pt")) +
  scale_colour_grey(start=0.2, end=0.8, name="Dependent variable") +
  ggtitle("")


explorative_plot_2

par(mfrow = c(1, 1))

explorative_plot_full <- ggarrange(explorative_plot_1, explorative_plot_2, ncol=1, align="v", common.legend=T) 
explorative_plot_full

#pdf("Explorative_plot.pdf")
#explorative_plot_full
#dev.off()

## Appendix

# Table A1 
m1_control <- lm.cluster(manipulationstjek ~ treat_corona + alder + kon + erfaring, cluster = 'kommune', data=deservingness)
summary(m1_control)

m2_control <- lm.cluster(corona_sanktion ~ treat_corona + alder + kon + erfaring, cluster = 'kommune', data=deservingness)
summary(m2_control)

m3_control <- lm.cluster(corona_hjaelp ~ treat_corona + alder + kon + erfaring, cluster = 'kommune', data=deservingness)
summary(m3_control)

# Table A2
summary(lm.cluster(corona_hjaelp ~ treat_corona + alder, cluster="kommune", data=deservingness))
summary(lm.cluster(corona_hjaelp ~ treat_corona*alder, cluster="kommune", data=deservingness))

summary(lm.cluster(corona_sanktion ~ treat_corona + alder, cluster="kommune", data=deservingness))
summary(lm.cluster(corona_sanktion ~ treat_corona*alder, cluster="kommune", data=deservingness))

# Table A3
summary(lm.cluster(corona_hjaelp ~ treat_corona + kon, cluster="kommune", data=deservingness))
summary(lm.cluster(corona_hjaelp ~ treat_corona*kon, cluster="kommune", data=deservingness))

summary(lm.cluster(corona_sanktion ~ treat_corona + kon, cluster="kommune", data=deservingness))
summary(lm.cluster(corona_sanktion ~ treat_corona*kon, cluster="kommune", data=deservingness))

# Table A4
summary(lm.cluster(corona_hjaelp ~ treat_corona + erfaring + alder + kon, cluster="kommune", data=deservingness))
summary(lm.cluster(corona_hjaelp ~ treat_corona*erfaring + alder + kon, cluster="kommune", data=deservingness))

summary(lm.cluster(corona_sanktion ~ treat_corona + erfaring + alder + kon, cluster="kommune", data=deservingness))
summary(lm.cluster(corona_sanktion ~ treat_corona*erfaring + alder + kon, cluster="kommune", data=deservingness))

# Table A5
summary(lm.cluster(corona_hjaelp ~ treat_corona + brugerorientering + erfaring + alder + kon, cluster="kommune", data=deservingness))
summary(lm.cluster(corona_hjaelp ~ treat_corona*brugerorientering + erfaring + alder + kon, cluster="kommune", data=deservingness))

summary(lm.cluster(corona_sanktion ~ treat_corona + brugerorientering + erfaring + alder + kon, cluster="kommune", data=deservingness))
summary(lm.cluster(corona_sanktion ~ treat_corona*brugerorientering + erfaring + alder + kon, cluster="kommune", data=deservingness))






